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ABSTRACT 


Numerical tables available for M/E,/c queueing systems 
are discussed. A new approximation method for steady-state 
information and waiting time distribution of this queueing 
system was developed. Validity of approximation was estab- 
lished directly for the large waiting times and by simula— 
tion for the smaller values. The developed method enables 
One to find delay probability, expected number in the 
queue and in the system, expected time to be spent in the 
queue and in the system, and probability of waiting for 


more than a specified time t. 
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I. INTRODUCTION AND SUMMARY 


Multi-server queues constitute a large proportion of 
queueing systems which arise in practice. Among those, 
queues with Poisson arrivals and Erlang service times (which 
W1ll be denoted as M/E,/C) occupy an important position, as 
Seen in the banks, airport check-in counters, hotels, super- 
Markets etc. However, no significant theoretical work was 
done until the late 60's. In recent years, methods for the 
Steady state solution of M/E,/c queues became available, fol- 
lowed by the numerical tables which were obtained by ee 
mentation of these methods. Even so, there is still some 
computational difficulty involved obtaining some particular 
information about above-mentioned queues, such as probabili- 
ty of delay exceeding a specified time length. 

The importance of the M/E,/C queue 1S due to the fact 
that it 1s used to model any queueing system whose service 
time distribution is believed to be unimodal. The solutions 
are known for extreme values of k: exponential service time 
distribution for k=1 and constant service times as k>, 
Usable solutions for multi-server queues having either expo— 
nential or constant service times are readily available. 

The M/E,./c queue is also important by providing a large 
variety of service time distributions in between these two 


extremes. 





The purpose of this study is, by means of computer simu- 
lation and examination of numerical tables published earlier, 
to present a simple and accurate computation method for ob- 
taining information about weestate solution of the M/E,/C 
queue. 

For this purpose, a general description of M/E,/c queue— 
ing model with definitions and assumptions for the analysis 
of it are given in Chapter II, followed by a summary of for- 
mer studies. 

In Chapter III, the numerical tables are analyzed and 
Some conclusions are drawn in terms of a simple form approxi- 
mation for steady—state distribution of waiting times. 

This approximation for the distribution of waiting times 
needed verification for small values of them. These cases 
were studied by computer Simulation. The procedure and the 
results of simulation are given in Chapter Iv. 

The results of the analysis of numerical tables and simu- 
lation are combined and generalized in Chapter V. Then the 
computational method developed by this generalization is de- 
scribed on a step-by-step basis. With this method, there is 
no need for numerical tables, computations are very simple, 
and the results are accurate for most practical purposes 
(the error being about 2% in general). To demonstrate the 
advantages of the proposed method, a comparison of it with 
the other method which uses numerical tables is given at the 


end of that chapter. 
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Pia rne M/E, /¢ Queue 


The description of the M/E, /C queueing model with assump— 
tions, definition of parameters and system variables along 
with some basic Bee ieee in the queueing theory forms the 
first half of this chapter. A brief discussion of the pre- 


vious studies on this model then follows. 


A. DESCRIPTION OF THE SYSTEM 

The M/E, /c queue is a multi-server queueing system with 
c servers, where arrivals occur according to a Poisson pro= 
cess with rate A, and service times have Erlang distribution 
With shape parameter k. It is also an infinite queue, i.e. 
there is no limitation for the number of customers in the 
system. 


The distribution of interarrival times X is given by 


EMeome emo nx > 0 


and the service times have the following probability density 


function: 


T ee 


Eety).= 
: (k-1) 19% 


with mean kg and variance ke". 


eat 





1. Definitions of Parameters 
a. Arrival Rate and Service Rate 
Arrival rate to the system is the reciprocal of 
mean interarrival time and denoted by \. Similarly, service 
rate of a server is defined as reciprocal of mean Service 


time and represented by pu. Then, u = = 


8 ® 
b. Traffic Intensity 


Traffic intensity p is defined by the following 


expression: 


0 = Afeu . 


c. Offered Load 
Offered load is the ratio of arrival rate A to 
service rate wu and denoted by a. It can also be expressed 


by product of number of servers c and traffic intensity o0 


as follows 


d. Coefficient of Variation 
The ratio of standard deviation to mean is de- 
PeWeduaseGOetricient Of Variation of a probability distribu— 


tion. It is denoted by V. For Erlang distribution coeffi- 


cient of variation is expressed as 
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2. System Variables 
a. Number in the System 
The total of the.number of customers waiting in 
the queue and those in the service is defined as number in 
the system and denoted by N. The number of customers only 
in the queue is represented by Nq. Expected values of N 
and Ng are denoted by L (average number in the system) and 
Lq (average number in the queue), respectively. 
b. Waiting Time 
The time spent in the queue by a customer is de- 
fined to be the waiting time of a customer. The time spent 
in ONew system is then the sum of waiting time and service 
time of a customer. Tq denotes waiting time and T denotes 
the time spent in the system. Expected values of T and Tq 
are denoted by W (average time in the system) and Wq (aver- 
age waiting time) respectively. 
c. State Probabilities 
The number in the system at a particular time 
is defined as the state of the system. Then state probabili- 
ty P(N = nit) is the probability of the system being in 
state n at time t. | 
d. Delay Probability 
Delay probability is the probability that the 
arriving customer finds all the servers busy and enters the 
waiting line. It is denoted by P(Tq > 0). Also, by defini- 


Eon Lt f£Ollews that 


ee 2 (Nez). 
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3. The Steady-State Solution 
The steady-state solution is defined to be the proba- 
bility distribution of the number in the system when the 
ete achieves statistical equilibrium. In the steady-state, 
the state probabilities do not depend ont, i.e. P(N = n|t) = 


P(N =n). Also it follows that 


oo 
» P(N=n)=l1. 
n=0 

In queueing theory, it is a well-known fact that [1] 
the steady-state is achieved if and only if the traffic inten- 
Sity p is less than l. 

4. Assumptions 

The M/E, /c queueing system under study is assumed to 
have the following properties: 

(ar) There is only one waiting line no matter how 
many servers there are. 

(11) The customer at the head of the waiting line 
starts getting service immediately whenever a server becomes 
available. 

(11i) Order of service is first-come first-served. 

(iv) If an arriving customer finds all servers busy, 
he joins the waiting line. The waiting line has unlimited 
capacity. 

(v) The service channels (servers) are indexed 


consecutively. If an arriving customer finds more than one 
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server vacant, he goes to the server with smallest index 
(this assumption doesn't cause any loss of generality, but 
is convenient for the computer Simulation model). 

(vi) The servers are homogeneous, i.e. the service 
distribution is the same for all servers. 

(vii) Only one customer at a time can be served by 


amoeCriver. 


Eee SORMER STUDIES 

There are numerous studies in the literature of queueing 
theory, done for the steady-state solution of M/E,/c queue. 
But most of them cover only some particular values of k and 
c, and their results cannot be generalized. Therefore, such 
Studies will not be mentioned here individually. 

The earliest suggestion known by the author about the 
solution of general M/E, /c System was mentioned by Lee [2]. 
He stated referring to a case study done in 1956 that mean 
walting time can be approximated by use of David Owen's sug- 
gestion. The formula given in [2] was said to be usable for 


tae VY < 1.0 and as follows, 


W = 5 W (1+v*) 


eal 


where : mean waiting time for M/E, /c queue 
: mean waiting time for M/M/c queue 


; coefficient of variation of Erlang (k) 


Geseribution, v2 = 1/k. 


The validity of the approximation was demonstrated by 


comparing its results with simulation results for different 


eS 


cases. However, the approximation is just for mean waiting 
time or queue length, hence it isn't possible to approximate 
state probabilities. Also, no reference was given for de—- 
tails of David Owen's suggestion. This approximation will 
be mentioned again later in this study. 

The steady—state solution of the general M/E, /c queue 
was first studied by Mayhugh and McCormick [3]. The results 
can be applied to the cases with any value of k and c. How= 
ever, the computation procedure iS so complex that it 
requires a very considerable amount of work. 

A parallel study was also done by Heffer [4]. The solu- 
tion method proposed by Heffer differs from Mayhugh and 
McCormick's, but it also is very complex. 

The reader is referred to the references [5] and [6] for 

more detailed discussion of these two studies. 

The most recent study was done by Yu [5]. It concerned 
the En! By/e queue with heterogeneous servers, but results 
were also stated for nomogeneous server case. Then letting 
m= 1 gives the solution procedure for M/E, /c queue. However, 
like the other two studies, the computations required still 
demand an enormous amount of work. 

The theoretical results obtained by Heffer [4] and Yu [5] 
formed the basis of the only numerical tables available for 
M/E,/c queue, prepared by Hillier and Lo [6]. The tables 
have cases of M/E, /c queue with limited values of k and c 
(Sec peCus Oye did ea Lew cases for Ey! Ey/e queue. These 


tables will be discussed in detail in the following chapter. 
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Ili. DISCUSSION OF NUMERICAL TABLES AND 
DEVELOPING A NEW COMPUTATION METHOD 
In this chapter, a detailed review and discussion of the 
numerical tables given in [6] will be made. Then, using the 
results of these discussions, a new computation method will 
be developed. 


Sample pages from the numerical tables are given in Appendix B. 


A. DISCUSSION OF NUMERICAL TABLES 
l. General Description 
The tables under consideration were designed for gen- 
eral El B/S systems. The analysis however will be focused 
On the cases in which m = 1 (i.e. M/E,/¢ system), which forms 
a large porportion of the tables. As mentioned earlier,. the 
values of k and c for various tables are as follows: 
Keteuc= lp 2, . age l0 
Korie tl Denes 5 
k=4; c=2,3 
Ka Opece pO; C=2 
The information given in each table is: 
(i) State probabilities, P(N=n), and 
(ii) Cumulative state probabilities, P(N<n), for 
ej usis.e 5 
(111) Delay probability, P(T,70); 


(iv) Expected number in the system, L; 


di] 





(v) Average queue length, La; 
(vi) Average queue length for the equivalent M/M/c 


a 


system with the same arrival and service rates, boi: 
(vii) The ratio L,/Lg1° 
All this information iS given in each case and for 
the following values of traffic intensity op : 
Pn OLMtOrmmOrrZ 0 sep Oe OO, On55,.-., 9.95, 0.98, 0.99. 
So, all the information is given for steady-state con- 
aeons . 


Average queue length L for the equivalent M/M/c sys- 


ql 
tem is computed by using the following formula [7], 


p= —el, p (3.1) 
= oR oy), 


where 





ea ; A 
it. > fep)e 4 _(ep) 
Po j=0 


The rest of the information given in the tables was 
computed by the results of theoretical work Mentioned in the 
last chapter, done by Heffer and Yu. 

2. Use of Tables 

The kind of information listed above can be obtained 
directly from the tables. However, by making use of some 
general relationships in the queueing theory, some other in- 
formation can be obtained besides those mentioned above. 

a. Interpolation 

No interpolation method is specified in the ex- 


planation sections for the tables in [6], for the values of 


dS 





0 different from those given above. The default position 
taken is linear interpolation. 
b. Waiting Times 
The average waiting time Wy and average time in 
the system W are frequently of interest in the analysis of 
queueing systems. These values can be obtained from the 
tabulated values of Je and L respectively, by using the fol= 


lowing relationships [8], 


ech = b 
Ng mie! W x 
The probability as that waiting time exceed— 
ing t is also important in the analysis of queueing systems. 
Although this kind of information is not given in the tables, 
P(T,>t) can be approximated from the state probabilities, 


P(N=n), as follows: 


eee = » Pa@i=n) PtP<[k(n-erl)—1) >; t>0, 
n=c : 


C2) 


where D is a Poisson random variable with mean ckut. The 
reader is referred to the discussion in [6] for the stochas— 
tic reasoning of this relationship. The quantity 

P{D< [k(n—c+1)—1]} can be obtained from Poisson distribution 


tables with mean ckut. If this mean exceeds 25, then the 


1s, 


normal distribution with the same mean and variance can be 
used as an approximation to Poisson distribution. 
3. Discussion | 
a. Delay Probabilities 

Similar to M/E,/C, M/M/c denotes a multi-server 
queueing system with Poisson arrivals of rate \ and exponen- 
tial service times of service rate uw. The delay probabili-~ 
ties for M/M/c queue can be computed by uSing the formula 


7, 


P(T >0) = _(cp)*_ - (nen 
q c!(l—o) ,c=l1 : is 
Ds (cp) 7 , {cp) 
4=0 3! c!(1—p) 


A chart is also available in Appendix A for c=2, 
Sewers and O<op<l, for P(T 79) - 

A comparison of delay probabilities of M/E,/c 
and M/M/c systems for selected values of o, ¢ and k is given 
in Table I of this study. A careful examination of this 
table shows that the corresponding delay probabilities of 
the two systems for the same value of p are very close. 
Since the formula or charts are available for M/M/c ayeten: 
tnis comparison shows that they can also be used to estimate 
delay probabilities for M/E,/C system with a very small 


error, usually less than 0.01. 
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b. Average Queue Length and Number in the System 
Pecateruleexamination of the "Ratio” column of 


Table I in [6] indicates that 
a al i 2 


which is essentially the same as Owen's suggestion [2]. How-— 
ever, it works beyond the limits of coefficient of variation 
V that were stated by Owen. Also (3.4) 1s true even when 0p 
is markedly less than one. The theoretical verification for 
(3.4) was developed by Yu in [9]. 

A‘'more precise approximation formula was also 
developed by Hillier and Lo in [6] which provides more accu 
racy for small values of po. In this case, the approximation 


formula is given by 


1 1 
L_ x x (1tg) (1 + i? L 


q (S45) 


gl 


where g = f(p,k,c). Then the expression for g has been ob—- 
tained by linear regression. The exact coefficients for 
g(po,k,c) are given in [6]. However, a rougher expression 
which is more convenient for computations will be used here, 


given as 


Ecal:, 
ae 


wv 


—1l 
Cran (c—1) 


magi aren (lo) -} 326) 


ZZ 





The main purpose of Hillier and Lo for introduc— 
ing g into approximation formula (3.5) was to provide means 
of approximation for the cases with larger k and c values 
which were not covered in [6], namely, extrapolation for lar- 
ger values of c and k. One way to check the validity of the 
Semputatlonal formula for g is to first let k go to infinity, 
then compare the values of average waiting time Wa computed 
by the formula given below with average waiting time ee Sis 
M/D/c queue. This latter queue has Poisson arrivals with 
rate X4 and constant service times of length 1/u (or kg), and 
One can use the fact that the distribution of service times 
can be approximated by constant service time distribution 
for very large k as mentioned in Chapter I. 


The computation formula for W of M/D/c queue 


qD 
is given in [2] as follows, 
co ©O 5 O° ; 
Z —ia (ia) 1 (1a) 
Wap ~ 2 : 2 “jr 7 5 2 Fr 
i=l jJ=1c over 


However, a more convenient form for computations can be ob- 


tained as 





(Sao 
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which permits the use of Poisson distribution tables. But 
it still requies a great deal of computational work. For- 
tunately, a very convenient chart was developed by Shelton 


eo] £or W for 1l<c<100 and 0.10<p<0.96. 


qD 


The comparison mentioned above can be made hav— 
ing the necessary tools available and Table II can be used 


for this purpose. The values of Wop for M/D/c system were 


obtained from Shelton's charts. Average waiting time NG for 


M/E, system 1s computed by using the following formula, 


273 1 


a, ok = ee 
w= ((1-p) + (1-9) 73] + Poy 


q [2 + + og lppp) (OL) 
where Lo 1s computed from (3.1). If the expression given 
by (3.6) is true for very large values of k, then it should 
also be true that 


lim W_.=W 
| gD 


k>0o 


Taking the limit of Wy gives 


| Ee 2/3 ae 
lim We = 31 + 1) (c—1) {(l—p) + (lp) 1] E 


a 


gl 


kK>0o 
On the other hand,if g 1s omitted, the limit 


will then be 


| i = 1 
at em Dy qlee sone” 


k>0o 


ql 
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gD and these two limits for 


eight different cases is given in Table II. Investigation 


A comparison of W 


TABLE II. 
Comparison of Wop and Limits of Wa With or Without g 
Case No Cc 0 = Moy Bg ee Wop 
1 2 0.80 tee e209 1.08 
2 2 0260 1.422 1.451 ee 6 
3 5 0.90 Zeere Ont 2.34 ZielO 
4 5 0.90 1.144 LL iOS 
5 10 0.90 er oeu7 Paco. i.36 
6 10 Oe 0 ie 0 5 Pe263 eS 
7 20 0290 3305 5 OG a6 
8 20 OS ENG 5.508 5.88 5.6 


of these results shows that, for almost all the cases the 


QD than the other limit. 


This indicates that the formula (3.6) should be reviewed 


limit with g omitted is closer to W 


for large k values. However, this revision can be done only 
after some additional tables become available for the cases 
not covered in [6]. 
c. Waiting Times 
The method of approximating the probability 
P(T,>t) that the waiting time will exceed t using the state 
probabilities P(N=n) was described earlier in the chapter. 


However, this method assumes that an arriving customer will 
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have to wait for approximately k(n-c+l) service phase (ex— 
ponentially distributed with mean 8) completions before a 
Server becomes available to start his service. It is stated 
in [6] that this approximation for p(T?) should be better 
for larger values of po and t due to this assumption. It 
should be added that, for large k and small t, the approxi~ 
Mation fails greatly in representing the actual distribution 


Of waiting times for the same reason. 


B. NEW COMPUTATION METHOD 

The construction of a practical computation method using 
the facts obtained from the preceding discussions will form 
the rest of the chapter. This new method is desired to be 
applicable to the queueing problems involving M/E,/C systems 
met in practice, with simple and quick computational work 
without any need for tables, interpolations etc. 

The kinds of information which are desired to be able to 
compute by the new method for M/E, /C queueing system are 
mainly 

(1) Delay probability, 

(11) Average waiting time, average time in the system, 
average queue length and average number in the system, 

(iii) Probability distribution of delay times and proba~ 
bility that waiting time will exceed a specified time length. 

1. Delay Probability 

In the discussion earlier it was pointed out that 


the delay probabilities P(T,>0) for the systems M/E, /c and 
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M/M/c are very close for the same value of p. An examina- 
tion of the comparisons given in Table I shows that the 
delay probabilities are very close or M/M/c (k=1), M/E, /c 
(k>1) and M/D/c (k=) systems with the same value of 0. 
The differences between the delay probabilities for the same 
o are less than 0.01 for almost all the cases and gets even 
smaller for p close to one. This comparison suggests that 
delay probabilities of the M/M/c system can be used for cor= 
responding M/E, /C system also. These delay probabilities 
are computed by uSing formula (3.3). The formula is not 
Suitable for large values of c, however, and a chart is given 
in Appendix A for up to about 16. For greater values of c, 
the charts are also available in [10] and [11]. 
2. Average Queue Length and Average Waiting Time 

In the previous discussion, after examination of the 

tables in [6] it was concluded that the average queue length 


for M/E,/c system can be approximated by 


No] 


i = 


Hl 
- eg el + RE gi (3.8) 


where el 1s computed according to (3.1), and g is computed 


by 


ys “> 


= 2/3 


g= & Sle)? {(1-p) + (1-p)*} . 


Contribution of g in equation (3.8) is negligible 
for practical purposes in most cases. Also it was shown 


earlier in this chapter that the formula g requires some 


Ze 





modifications. Therefore it will be omitted in further 
discussions. 

Now, attention will be focused on equations (3.1) 
and (3.3). Rewriting equation (3.1) ina different form, 


then substituting (3.3) gives 


L = Op. (co) ~ and 
Saare > te)) . tee) 
a cp cp 
c! (l=) ad ie + Bir Cies 
o°P(T _>0) 
if 


ql (1-9) 


Then it follows that 


Lai P(T >0) 
W = = 
ql d oy) (fe) 
Using (3.8) with g=0, some approximation formulas 
for average queue length and average waiting time are ob- 


tained respectively as follows 


P(T >0Q) 


_1 1 

L, = 50 qa lt (3.9) 
1 P(T _>Q) l 

aecemcnaaor °f eer 


® 


where P(T,79) 1s delay probability for M/M/c system which 
can be obtained from the charts very easily. Then it is 
very simple to compute Ly and Wa uSing equations (3.9) and 


oie 0). 
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The average number in the system L and average time 
in the system W can be computed by following well—known rela— 


tionships of queueing theory, 





3% Waiting Time Distribution 


Suppose that the conditional waiting time distribu- 
tion given 1,79 ON M/E, /c queueing system can be approxi- 


mated by 


p(T >t|T>0) =e >° , 
q q 


Then uSing Bayes' theorem it follows that 


—bt 
P(T >t) =P Se}. SO > a Oy % Bel 
( q ) (T, IT DENTS 0) e ac ) ( 1) 


Let Pog be the CDF of waiting times. Then 


Pree. —bt 
Ppg(t) = 1-P(T,>0) e ; t 0. 


Now average waiting time can be computed as follows 


_ mer ke _— —bt 
Haar See (ime (t}) dt = s P(T >0)e vat 
P(T >0) 
Wo = —i— : (3.12) 


Suppose the parameter b 1s modeled as 


= 2cu(l—p) 
(1+V*) 
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Then, 
P(T >0O) 


_ 2 
q  2cu(1-p) oe) eae 


W 


where V is coefficient of variation of service time distri- 
bution. Notice that equation (3.13) is exactly the same as 
equation (3.10). However, there can exist some other wait— 
ing time distribution to give the same average waiting time 
as given by (3.13), so it must be shown that for M/E, /c sys— 


tem, the distribution of waiting times can be approximated 


by 


= 24coia eg hs 
Fog(t) = 1-P(T>O)e ~(1+V2y——i«‘««#;:«é (3.14) 


Pe (3el4) 25 true, then ee es computed by (3.11) 
should be approximately equal to the one obtained by the pro— 
cedure suggested by Hillier and Lo [6] which uses state pro- 
babilities as described earlier by equation (3.2). Table 
III gives the comparison of ee values obtained by these 
two different formulas for various values of t. It should 
be recalled that approximation by state probabilities is 
poor for small values of t. However, for caw the two re- 
sults agree very well. The difference is ease of computa- 
tions. Equation (3.11) is much easier than equation (3.2) .- 
to compute the same value. 

The validity of approximation (3.14) for small values 
of t will be shown by simulation, which forms the next chap- 
ter. Then the suggested computation method will formally 


be described on a step-by-step basis in last chapter. 
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IV. COMPUTER SIMULATION 


A computer simulation model for M/E,/c System is used 
to investigate the validity of approximation for waiting 
time distribution introduced in the preceding chapter. 
Validity of this approximation for moderate and large 
values of t was demonstrated in the same chapter. There— 
fore, analysis will now be focused on the approximation 


for small values of t. 


A. COMPUTER PROGRAM 
1. Model 

Since the model used to construct the simulation 
program is based on the same assumptions as stated in Chap-— 
ter II, it won't be described again in this chapter. The 
computer language selected was SIMSCRIPT which is especial- 
ly convenient for simulation of multi-Server queueing sys- 
tems. The reader is referred to [12] and [13] for a brief 
idea about SIMSCRIPT if he is not familiar with this 
language. Main reference for SIMSCRIPT however is [14]. 

One way to check for small t the validity of approxi- 
Mation distribution developed earlier is to keep the fre- 
quency table of waiting times of the customers during the 
Simulation. Let M(t) be the number of customers with wait- 


ing time greater than t and M be total number of customers 
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> 


served during simulation period. Then the estimator for 


probability of waiting time greater than t is 


—.., 


P(T,>t) = M(t) 


M 





ae.) 


This value can be compared with the value given by the ap- 
proximation formula to check its validity. 

The computer program is given in Appendix B. It is 
written in the way that it would be possible to obtain the 
estimates of all kinds of information about the queueing 
system being simulated; average waiting time and average 
gueue length, state probabilities, delay probability etc. 
However, only the estimates given by (4.1) will be discussed 
here. 

2. Variance Reduction 

One of the problems encountered in the simulation 
of queueing systems is the high variability of waiting times, 
especially when traffic intensity p is high. Also strong 
positive correlation between the waiting times of consecu- 
tive customers causes serious errors if the standard for- 
mula is used to estimate the sample variance of waiting 
times since this formula will underestimate the true vari- 
amee i). 

There are several methods developed to overcome this 
difficulty with high variability. Interested reader can find 
comprehensive information in references [12], [15] and [16]. 
The variance reduction method which will be used here is 


antithetic variates. 
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In the simulation program, two different random num— 
ber streams are used to generate interarrival and service 
times. Once a uniform variate between 0 and 1 is generated, 
then the random variate from a desired distribution can be 


Obtained by 


X =F (u) (Blegy) 


where u is a uniform variate between 0 and 1 and aia denotes 
the inverse of CDF of X [12]. Then, for example, a sequence 
Of exponentially distributed random variates with mean 1/i 


can be obtained from a stream of uniform variates by using 


am 


X = ln (u) (4.3) 


Let Zy and Zo be estimators of a parameter, obtained 
from two different Simulation samples and Z3 = x (2, + Zo) 


Then 


+ 


- aia): 
Var (Z.) = gq Var(2Z,+2Z.) =r [Var (Z,)+Var(2,)+Cov(Z,,2Z,) ] 


(4.4) 
which implies that maximum negative correlation between Zy 
and Zo would minimize Var (Z3). This can be obtained by us- 
ing l—u in (4.2) in place of u for random variate genera- 
tion in the second simulation run. 


During the simulation, a random sequence of large 


service times or short interarrival times would cause long 
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waiting times and conversely. Using the procedure described 
above, a sequence of large Peeing times would become sequence 
of short waiting times in the second run; in other words, the 
Baiting times in two different runs will be negatively corre— 
lated. Let P, (t) be the probability that waiting times will 
exceed t in the ith run, 1= 1,2. Using antithetic variates, 
P, (t) and P(t) will be negatively correlated; then the esti- 
mator will be computed by 


Sw 1 
PATS t) = (Py (t) + P.(t)) 


where P, (t) and Po (t) are computed according to Aan 

Generation of antithetic exponential variates for in- 
terarrival times is given’by (4.3). However, it 1S not 
straightforward to generate antithetic Erlang variates since 
the CDF cannot be inverted in closed form. Generating Erlang 
variates as sum of k exponential variates with mean 8 does 
not help since equation (4.2) cannot be utilized to obtain 
antithetic variates. Nevertheless, one way to solve this 
problem is to store the CDF table of chi-square distribution 
in an array, then compute Fo *(u) by linear interpolation to 
get a chi~square random variate, and obtain the Erlang 


variate by the following transformation 


where E, 1s the Erland variate with mean k&g and Cy is the 
chi-square variate with degrees of freedom k. The tables 


for CDF of chi-square distribution are available in [17]. 
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Be yWace Collection 
a. Selection of Input Parameters 

The input values for each pair of runs were arbi- 
trarily selected to give certain traffic intensity p values, 
usually 0.8 and 0.9 in order to get a wider range of waiting 
times even though they cause high variability. The width of 
frequency intervals was chosen to be 0.25 min. so that fie 
for example, would show the number of customers with waiting 
time less than 0.25 minutes. Then M(t), the number of 


customers with waiting time greater than t can be computed 


by 


co 


M(t) = »s f 


1=4t+1 


Since t will be in multiples of 0.25. This computation is 
done in the computer program. 
b. Initial and Final Conditions of the System 

Initial and final conditions of the simulated 
system are important since they effect the value of para> 
meter estimated. The approach taken in this study was to 
use epochs. If the time at which the number in the system 
N changes from 0 to 1 by an arrival is defined as a regenera- 
tion point, then the interval between two successive regene— 
ration points can be defined as an epoch [12]. It is 
Obvious that the sequences of waiting times in two different 
epochs will be independent from each other. The epochs tend 


to be lengthy with high traffic intensity o, large arrival 
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rate A, with large number of Servers c, or a combination of 
these three factors. Before each simulation run, the number 
of epochs for which the Simulation is to be run was deter— 
Mined as an input conSidering those three factors. Exper- 
lence showed that the first few epochs tend to have waiting 
times smaller than usual, so they were omitted for data 
analysis. This isn't feasible for cases with long epochs; 
however the first couple of epochs have enough information. 
The simulation runs were started with N = O and an arrival 
so that starting time would be a regeneration point and 
Stopped after the number of epochs determined earlier is 
completed, i.e. N = 0 was the final condition. 

P, (t) or Po (t) are given in the output of program 
for t = 0.25, 0.50, 0.75, .... etc. Then to get the estima- 
tors P(T,7t) all one has to do is to average them for cor- 


responding t values. 


B. RESULTS 

A comparison of waiting time probabilities obtained from 
Simulation and approximation method is given in Table IV. 
Investigation of this table indicates that the approximation 
method agrees quite well with the results of simulation for 
Small t values. The difference between the corresponding 
values for moderate or large values of t for which the 
validity of approximation was demonstrated earlier explains 


the difference between the values when t is small. 
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Ve. CONCLUSION 


It was shown in the last chapter that the approximation 
method developed earlier for the distribution of waiting 
times can also be used for small values of t. The method 
for computing steady-state information of M/E,/c queueing 


System can now be formally described. 


A. DESCRIPTION OF THEPROPOSED COMPUTATION METHOD 

After analysis of data and the decision has been made 
that the particular system under study can be approximated 
by M/E,/c system as described in Chapter II, and also having 
the estimators for i, k and 8 obtained (by maximum likeli- 
hood or method of moments), one is ready for computations 
to get the desired information for the system. 


1. Delay Probability 


First compute service rate u as 


ace. 
Hu kB “ 


Then compute traffic intensity p and offered load a by 


0 = — and a = 


e|> 


respectively. Then enter the chart in Appendix A with a to 
get delay probability P (7,79) - Use the charts in [10] or 


{11] if number of servers c>l6. 
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23 Average Waiting Time and Average Queue Length 


First compute a system constant, namely b, by 


Zewtl=p) 


L 
L+¢F 


b= 


Then average waiting time Wg 1s computed as follows 


P(T _>0) 
ee 
g b 


Also, average queue length Ly 1s given by 


L = *AW_. 
q q 


3. Average Time in System and Average Number in System 


Average time in system W and average number in sys- 


tem L will be computed by 


W= W_ + i and L= ZL +A 
Hy H 


respectively. 
4. Waiting Time Distribution 


The CDF of waiting times is given by 


tae —bt 
Fpg(t) = 1-P(T poe, t>0 


where b is the system constant computed in step 2. 


Then the following probabilities can be computed 


(i) P(T >ty) = P(T >0)e*1 

e q . q —bt -bt 

Gite <ta)i = p(T >0) (e l1l-e 2) 
age g —bt 


(25101: ) P(T<t)) 1-P (T,>0)e nu 


as needed by the user. 
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B. DISCUSSION 
i advantages Of the Method 

(1) Needless to say, the most important advantage 
of the method is simplicity. The kind of information listed 
above can be computed in minutes given c, A, k and uw which 
would be needed to compute anyway. All of the computations 
Can be laid down on one regular size page so that it is 
very easy to follow by somebody else. 

(ii) In most cases, one has to compute the above 
listed information to determine the effect of changing the 
number of servers c on the selected system variables. This 
multiplies the savings of time and computational effort. 

(iii) No tables are necessary except for the delay 
probability chart. No interpolation would be needed. 

2. Disadvantages of the Method 

(1) Since the method gives an approximation, it 
1s not too precise even though the results are almost always 
in + 5 percent of the actual value, and in + 2 percent for 
most cases. 

(11) The approximation of state probabilities with 
this method does not appear to be direct. 

3. Applications 

One of the characteristic applications of the queue— 
ing theory is to investigate the effect of changing the 
parameters or number of servers on the measure of effective- 
ness (MOE) assigned by the management. This measure of 


effectiveness can be delay probability, average waiting time 
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or probability of waiting time exceeding time t etc. The 
marginal utility of c + 1st server will be the difference 
between the values of MOE's for c and c+l servers. Then the 
decision criterion whether or not to hire the c+lst server 
is how his cost compares to his marginal utility. 

Another managerial objective may be to require, for 
example, the probability of waiting time exceeds 4 min., 
ae) = 0.05. The number of servers necessary to achieve 
this purpose will then be of interest. 

It is possible to extend the variety of examples. 
The common property of the above examples is that precision 
greater than that of suggested approximation is not needed 
to make the decision. 

4, Final Remarks 

The author wishes to emphasize that without the 
numerical tables provided by Hillier and Lo [6], it would 
have been impossible to develop such an approximation method. 
The method may need some modifications as new tables for 


the cases not covered in [6] become available. 








APPENDIX A 
CHARTS FOR DELAY PROBABILITY 


In the two following charts, the delay probabilities 
can be found for M/M/c system which were shown to be very 
close to those of M/E, /c system, with number of servers c 
up to about 16. One has to enter to chart with offered 
load a as abscissa, then delay probability P(T,79) is the 
ordinate value of the point on the proper c curve which 


has a@ as abscissa value. 
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